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Abstract. - We study numerically the nonequilibrium dynamical behavior of an Ising model with 
mixed two-spin and four-spin interactions after a sudden quench from the high-temperature phase 
to the first-order phase transition point. The autocorrelation function is shown to approach its 
limiting value, given by the magnetization in the ordered phase at the transition point, rric, through 
a stretched exponential decay. On the other hand relaxation of the magnetization starting with 
an uncorrelated initial state with magnetization, rui, approaches either rric, for rm > 0.5, or zero, 
for rrii < 0.5. For small nii and for rUi slightly larger than 0.5 the relaxation of the magnetization 
shows an asymptotic power-law time dependence, thus from a nonequilibrium point of view the 
transition is continuous. 
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Introduction. — Most of our knowledge about the 
properties of phase transitions has accumulated for con- 
tinuous or second-order transitions for which the concepts 
of scaling and universality, as weU as the application of 
the method of renormahzation group, have provided a 
deep understanding [1]. Considerably less is known about 
singularities at first-order phase transitions, although in 
nature this type of transition is very common [2]. At a 
first-order phase transition point one observes phase co- 
existence where several thermodynamical quantities, such 
as the internal energy and the order parameter (magne- 
tization), have a discontinuity, whereas the correlation 
length is generally finite. In spite of this finite correla- 
tion length some response functions, such as the magnetic 
susceptibility, have a scaling behavior in a finite system [3] 
that involves the discontinuity fixed-point exponent yd = d 
where d denotes the dimension of the system [4] . From a 
dynamical point of view the relaxation time is finite at 
a first-order transition point, yielding equilibrium auto- 
correlations and relaxation functions that decay exponen- 
tially. 

Many studies have been devoted to the investigation of 
the formation of the ordered phase when the temperature 
is lowered below the phase transition point, Tc- Thus, 



the nonequilibrium dynamics of the nucleation process [5] 
and the morphology of the solidification [6] are key prob- 
lems in material science. From a more theoretical perspec- 
tive one is interested in the coarsening process [7] that 
takes place when the system is suddenly quenched from 
the high-temperature phase at or below the phase transi- 
tion point. If the quench is performed below then the 
competition between the stable ordered phases leads to 
the same type of qualitative picture independently of the 
order of the transition at Tc [8]. If, however, we quench 
the system at the very transition point, important differ- 
ences are expected in the two cases due to the different 
behavior of the correlation length. 

In the past this type of phenomena, i.e. nonequilibrium 
critical relaxation, has been thoroughly studied in the case 
of a second-order transition point [9,10], whose investiga- 
tions were also motivated by the appearence of ageing. On 
the contrary only little attention has been paid to the case 
when the transition is of first order. It has been demon- 
strated that nonequilibrium relaxation starting from an 
ordered or a mixed-phase initial state can be used to lo- 
cate accurately the phase transition point and to decide 
about the order of the transition [11,12]. For models with 
quenched disorder the tricritical value of the dilution has 
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also been studied by this method [13]. However, we are 
not aware of any study in which the singularities in the dy- 
namical quantities are systematically investigated. In this 
paper we perform an investigation of the nonequilibrium 
dynamics at a first-order phase transition that addresses 
questions well established for a second-order transition 
point. In particular we consider the asymptotic behavior 
of the autocorrelation function as well as the relaxation of 
the order parameter after a sudden quench from the high- 
temperature phase to the first-order transition point. 

The structure of the paper is the following. After a 
brief recapitulation of the known results of nonequilibrium 
critical relaxation we present a detailed numerical study 
performed on a square lattice Ising model with mixed two- 
spin and four-spin product interactions [14,15]. This sys- 
tem has a first-order phase transition and the location of 
the transition point is known by duality. The numerical 
results obtained in this model are consistent with the ex- 
istence of a divergent nonequilibrium relaxation time. We 
interpret its origin and argue that the observed nonequi- 
librium dynamical behavior is generic at first-order phase 
transitions. 

Nonequilibrium critical dynamics: a reminder. — 

At a second-order transition point both the correlation 
length, ^, and the relaxation time, r, are divergent, and 
they are related through: r ~ The dynamical expo- 
nent, z, which generally depends on the local dynamics, on 
conservation laws, and on symmetries, is enough to clas- 
sify the dynamical universality class at equilibrium [16]. In 
a nonequilibrium situation, when the system is quenched 
from the high-temperature phase as an initial state to the 
critical point, new nonequilibrium exponents have to be 
defined [17, 18]. The reason for this is the broken time 
translation invariance due to a discontinuity at the time 
horizon ("time-surface"). As a result the autocorrelation 
function. C(s,t). is generally non-stationary: it depends 
on both the waiting time, s, and the observation time, t. In 
the limit t » s we have asymptotically: C{s,t) ~ 
which involves the nonequilibrium autocorrelation expo- 
nent A [18]. In another related process the initial state is 
prepared with a small, non- vanishing magnetization, nii, 
and one measures its relaxation, which for short times be- 
haves as m{t) ~ rriit^. Here the initial slip exponent 9 
satisfies the scaling relation: X = d — 6z [9, 17]. We note 
that in mean- field theory 9 = 0, whereas in real systems 
the exponent 9 is generally larger than zero and can be 
expressed as ^ = (xt — x)/z where x = fi/u is the anoma- 
lous dimension of the bulk magnetization (/? and v are the 
standard magnetization and correlation length exponents, 
respectively) and Xi is called the anomalous dimension of 
the initial magnetization. 

Ising model with multispin interaction. At a 

first-order transition point the correlation length stays fi- 
nite, ^ < oo, and the same is true for the equilibrium 
relaxation time, r < oo. The magnetization displays a 
jump from zero above to some value rric > below 



Tc. In this case the static critical exponents are formally 
given by x = p = and i/ = 1/rf, where the latter is the 
discontinuity fixed point value [4]. The specific model we 
study numerically in the following belongs to a class of 
square lattice Ising models with two-spin interactions in 
the vertical direction and n-spin product interactions in 
the horizontal direction that is defined by the Hamilto- 
nian 

-f^ = ~ X/ I -^ao-ija-jj+i -I- J„ CTj+fcj I . (1) 

The model is self-dual [14,15,19] and the self-dual point is 
located at: sinh(2J2/fcTc) sinh(2J„/fcTc) = 1. Forn = 2 
we of course recover the standard Ising model. Whereas 
for n = 3 the transition is found to be continuous and to 
belong to the universality class of the four state Potts 
model [20-25], for n > 4 the transition is of first or- 
der [26,27]. In the following we consider the model with 
four-spin interaction, in which case the ordered phase has 
eightfold degeneracy [28] and the phase transition of the 
system is comparable with that in the eight-state Potts 
model. If one takes J2 = J4, numerical results [27] in- 
dicate a latent heat of A/kTc = 0.146(3) and a jump of 
the magnetization from zero to nic = 0.769(6). The cor- 
relation length has not been measured in this system, but 
from the snapshots of the simulations one expects some- 
what different correlation lengths in the horizontal and in 
the vertical directions that are both of the order of ten 
lattice constants. 

Single spin autocorrelation function. The sin- 
gle spin autocorrelation function of the model is defined 
as C{s,t) = {aij{s)aij{t)), which is then averaged over 
all spins, 1 < i. j < L, and periodic boimdary conditions 
are applied in both directions. Starting with a disordered 
initial state, which corresponds to T = 00, we apply non- 
conserving spin- flip dynamics at the temperature T = Tc- 
In the following we consider no waiting time, i.e. s = 0, 
and use the notation: C{0,t) = C{t). In order to moni- 
tor finite-size effects we simulated systems containing 
spins with L ranging from 60 to 640. For the largest size 
we averaged over 2000 different samples with different re- 
alizations of the noise, more samples were simulated for 
smaller system sizes. As usual, time is measured in Monte 
Carlo steps where on average every spin is flipped once 
per step. The calculated autocorrelation functions for dif- 
ferent finite systems are plotted in Fig. 1. Rather large 
systems are needed in the calculation in order to han- 
dle the finite-size effects. From Fig. la we can identify 
three different time-regimes, i) For short times, which are 
shorter than the equilibrium relaxation time, t < t, there 
is a fast decay of C{t), which is described by a stretched 
exponential function, ii) This initial decay is followed by a 
relatively slow approach to the equilibrium limiting value 
limt^oo limL_+oc C'(i) = fn-c- Indeed the plateau of C{t) 
in Fig. la gives an accurate estimate nic = 0.770(3) in 
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Fig. 1: (a) The autocorrelation as a function of time for systems composed of spins with L — 60, 90, 120, 240, and 640. Three 
different time regimes can be identified, see main text, (b) The connected autocorrelation function C{t) — rric as a function of 
time. A stretched exponential decay (dashed line) is observed for times t > t with r ~ 200. 



accordance with previous numerical results [27]. For the 
largest finite system we have checked the behavior of the 
connected autocorrelation function and we have obtained 
a stretched exponential form: C(t) — mc ^ exp(— 6^"), 
with an exponent a — 0.198(3), see Fig. lb. We note 
that in the short-time regime the decay is described by a 
somewhat larger effective exponent in the stretched expo- 
nential, d « 0.32. iii) The third regime of C{t) is seen in 
finite systems [29] , in which it starts to decrease from the 
saturation value after passing a cross-over time: ^ i^. 
The dynamical exponent measured in this part of the fig- 
ure is found to be z = 2.05(10), which is consistent with 
z = 2. 

At this point we can conclude that for large systems the 
nonequilibrium autocorrelation function at a first-order 
transition point approaches the bulk magnetization, m^. 
Thus although the average magnetization of the system 
is vanishing (see the next section) the dynamical correla- 
tions display a finite asymptotic value. On the other hand 
the connected autocorrelation function has an asymptotic 
stretched exponential decay. 

Relaxation of the magnetization. Next we in- 
vestigate the relaxation of the magnetization at Tc where 
we start with an uncorrelated initial state with a magne- 
tization rrii < 1. In order to obtain reliable data for this 
noisy quantity, we averaged over typically 250000 different 
samples with different initial states and different realiza- 
tions of the noise. The data discussed in the following 
have been obtained for systems containing 240^ spins. We 
checked that the presented data are free of finite-size ef- 
fects by simulating some larger systems for selected values 
of the initial magnetization m^. 



The results are shown in Fig. 2 for different values of 
the initial magnetization. Here one can observe two fixed 
points concerning the limiting value of limt^oo rn(t) = 
moo, see Fig. 2a and 2b. For large enough initial order, 
nii > TO*, the magnetization approaches its finite equilib- 
rium value in the ordered phase, TOqo = TOc, whereas for 
weak initial order, rrii < to*, the limiting magnetization 
is vanishing: uioc = 0. According to our numerical es- 
timates the border between the two regimes is given by 
TO* = 0.500(1), and starting the relaxation from this ini- 
tial value the magnetization approaches a finite limiting 
value: TOqo = 0.3844(3), which is compatible with mc/2. 
At this special value of rrii there is first a drop of the 
magnetization in the early times when nuclei of the co- 
existing phases are created (resulting in an average mag- 
netization of « TOc/2), which afterwards takes part in a 
special coarsening process during which the coexisting or- 
dered and disordered phases are separated and thus the 
interface between the two phases are reduced. The eval- 
uation of the structure of the system in time is shown in 
Fig. 3. Here we present also results for to^ <^ in*, when the 
mass of the ordered domains and thus the magnetization 
in the system is continuously decreasing. 

In the regime rrii > the t dependence of m{t) is 
non monotonic when rrii is close to to*: in the early 
time-steps there is a sudden drop of the magnetization 
due to the formation of disordered clusters, which is then 
followed by an increase of m{t) that results from the 
formation of ordered domains. In this regime the re- 
laxation for t > T is described in a power-law form: 
m{t) — mc/2 ~ (rrii — m*){t — r)^. Here the measured 
effective exponent has a strong rrii dependence and close 
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Fig. 2: (a) Time evolution of the magnetization when starting from a state with an initial magnetization rrii, with rrii = 1, 
0.8, 0.7, 0.6, 0.51, 0.5, 0.49, 0.4, and 0.2 (from top to bottom). The grey line separates the two different fixed point, (b) The 
same as (a), but now the magnetization is shown as a function of Int, thus illustrating the non-monotonic behavior that one 
may encounter at early times, (c) m(t)/mi as a function of \nt for different small values of rrii. The initial decay is given by 
0.7 — O.llnt. The slope of the dashed line is 0.1. (d) At later times the decay for small initial magnetization is given by a 
power-law with an exponent 9' ~ —1 (dashed line). Shown are data obtained for rrii — 0.2, 0.1, and 0.05 (from top to bottom). 





Fig. 3: Snapshots of the system starting with an uncorrelated initial state having an initial magnetization rrii ~ 0.5 (upper 
figures) and rrn =0.1 (lower figures) at times t — 5, t — 100 and t — 500, from left to the right. The four-spin interaction is 
given in the vertical direction. 
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to rui = m* it seems to approach 9 = 0.50(3). Start- 
ing with the maximally ordered initial state with nii = 1 
the relaxation of the magnetization is monotonous and fol- 
lows very closely the decay of the autocorrelation function, 
C'(t), which is measured by preparing the system in the 
fully disordered initial state, see Fig. 1. As for the au- 
tocorrelation function the approach of the limiting value, 
'nioo = Trie, is asymptotically given in a stretched expo- 
nential form: m{t) — nic exp{—b't°- ), and the measured 
exponent is a' = 0.196(5), which coincides with a within 
the error of the calculation. Wc can thus conclude that 
the exponents for the two functions are probably identical 
and they are given by: a = a' = 0.197(4). Note that at 
a second-order transition point {nic = 0) the decay with 
m, = 1 is given by m{t) ~ t~^^^, and our result does 
not fit with a naive application of the discontinuity fixed- 
pont value, a; = 0. In the regime with rrii < m* and for 
rrii close to m* the decay seems to be symmetric to the 
relaxation process observed for rrii > m*- On the other 
hand for a small rrii/m* < 0.2 the data points scales with 
m{t)/mi, and here we can observe also two different relax- 
ation regions. As shown in Fig. 2c, for short times, t < t, 
the decay is very slow, which is well fitted by a logarith- 
mic dependence: m(t)/mi ~ c — dint, with c = 0.7 and 
= 0.1. In the long-time limit, t > t, the decay becomes 
faster and is described by a power-law form: m{t) ~ , 
with 0' Ri —1, as shown in Fig. 2d. 

Discussion. — Nonequilibrium relaxation at a first- 
order transition point has considerable differences with 
the case when the quench is performed to a second-order 
transition point. This is due to the facts that i) the cor- 
relation length is finite, which leads to a finite initial time 
regime, t < t, and more importantly ii) in the equilibrium 
state there is phase coexistence between an ordered phase, 
having rric > 0, and a disordered phase with vanishing 
magnetization. In the relaxation process the system typi- 
cally evolves towards one of these phases depending on the 
magnetization of the uncorrelated initial state, rUi. Just 
at rrii = 0.5, which separates the two regimes, the system 
evolves into phase coexistence. At this special value of rrii 
the ordered and disordered phases are separated during 
a special coarsening process, thus the average magneti- 
zation approaches rrioo = rnc/2. If the initial magneti- 
zation is slightly above the border-line value, rrii = 0-5, 
after the initial drop of the magnetization the ordered 
phase starts to grow aroimd the interface of the extra 
droplets of the ordered phase. Here we use formally the 
scaling relation valid at a second-order transition point, 
= (xi — x)/z, in which we identify the anomalous di- 
mension of the initial magnetization with the dimension 
of the interface, Xi = d — 1, whereas we set a; = at 
a first-order transition. Consequently we have an initial 
slip exponent: = [d — l)/z, which is compatible with 
the numerical results in = 2. In the other case, when 
the initial magnetization is small, then in the initial state 
there are no long-range correlations, thus the system has 



a mean-field character. Since the transition point of the 

mean- field model is Tc{mf) > Tc, an effective coarsening 
process takes place during which the magnetic moments 
of the domains, n, grow practically with their volume, 
^ ^ rrii K^Y ^ f^^^ ■ As a consequence the total magneti- 
zation of the system decays very slowly with time, slower 
than any power of time. This coarsening goes on until the 
linear size of the domains l{t) ^ t^l^ approaches the corre- 
lation length ^. For longer times the magnetic moment of 
the domains stay constant, \x* ^ rrii, so that the magneti- 
zation start to decay as mit) ^ ir'^l'^ , which is compatible 
with the measured value in d = 2. 

The observed power-law dependences of the magnetiza- 
tion are somewhat unexpected, since it means that there is 
no relevant time scale in the problem. Thus in this sense 
we have critical nonequilibrium relaxation, even though 
the equilibrium transition is of first order. To understand 
this phenomenon we draw an analogy between nonequilib- 
rium relaxation, i.e. the time dependence of the magne- 
tization, m(t\ after a quench from T = cx) at i = 0, and 
the behavior of the equilibrium magnetization, m{rj), in a 
semi-infinite system measured at a distance, y, from the 
free surface located at y = 0. In the two systems time, 
t, and distance, y, play analogous roles and they can be 
related o&t ^ rf- , and similarly: t ^ . In both systems 
translational invariancc is broken at t = and y = 0, 
respectively. At a second-order transition point the equiv- 
alent initial-slip behaviors arc given by: m.{t) ^ fi^i-^)/^-^ 
and m(y) ^ yi^"^^)^ respectively, thus the anomalous di- 
mension of the initial magnetization, Xi, and the anoma- 
lous dimension of the surface magnetization, Xg, are anal- 
ogous quantities [30-34]. There is also an analogy when 
starting with a fully ordered initial state, m{t = 0) = 1, 
and having an ordcax^d siirfa('(\ rn{y = 0) = 1, when 
the asymptotic decays arc given by m(t) ^ fr'^l'- and 
"^(?/) ^ 1 respectively. Now let us see how these analo- 
gies work if the transition in the system is of first order. 
The properties of the surface transition in systems having 
a first-order transition in the bulk have been studied in the 
literature [35-40] . Interestingly, the surface transition is of 
second-order for weak surface fields, hs < hi, which turns 
to first-order for hs > h^, and the two regimes are sepa- 
rated by a tricritical point. This phenomena for hg < is 
known as surface induced disorder and is analogous to the 
wetting transition [41]. Indeed in the surface region the 
correlation length is divergent, and one should even define 
two diverging correlation lengths, and which are re- 
lated as 

^ with an anisotropy exponent, ^ > 1. Now 

turning back to the nonequilibrium relaxation process wc 
can interpret our results for small initial magnetization 
and for |mi — 1/2| <C 1/2 as due to a quench induced 
disordering effect, so that the nonequilibrium dynamics 
is critical, although the equilibrium phase transition is of 
first order. A closer analogy with surface induced disorder 
is seen if the relaxation starts with rrii = 0.5, where the 
system, after the initial short time decay, remains at the 
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magnetization TOoo < rric for any finite times. We note 
that this phenomenon is different to other kinetic analogs 
of surface induced disorder [42] and that this type of anal- 
ogy does not apply to the case of an ordered initial state 
and an ordered surface. The static magnetization profile 
at the ordered surface shows an exponential decay to TOc, 
whereas in the nonequilibrium relaxation problem the de- 
cay is found to be stretched exponential. 

The results obtained in this paper about the specific 
model in Eq.(l) are presumably generally valid for other 
systems having a first-order phase transition. In partic- 
ular the functional form of the relaxation function close 
to mi = and rrii = 0.5 are assumed to be described by 
universal exponents, which could depend only on the di- 
mension of the system. In this respect the type of dynam- 
ics (conserving and non-conserving) does matter, e.g. for 
conserving dynamics the dynamical exponent should be- 
come z = 3 instead of 2. On the other hand the stretched 
exponential form of the autocorrelation function is prob- 
ably system dependent. It would be interesting to study 
these questions also for other models. 
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